DECI + PREP + FAMI + ORAL + WRIT + PHYS,
data = USJudgeRatings))
summ(reg2)
?export_summs
export_summs(reg2)
# In order to export your regression tables in a clean format
export_summs(reg1, reg2) # to have publishable Table output
# export to doc or pdf or Excel or HTML
export_summs(reg1, reg2, to.file = "docx", file.name = "test1.docx")
export_summs(reg1, reg2, to.file = "xlsx", file.name = "test2.xlsx")
data <- USJudgeRatings
y <- data[,12]
x <- as.matrix(data[, -12])
export_summs(reg1, reg2, scale = TRUE,
error_format ="[{conf.low}, {conf.high}]") #play with the options
summ(reg2)
data <- USJudgeRatings
y <- data[,12]
x <- as.matrix(data[, -12])
###### if you want to remove more than one variable:
x2 <- as.matrix(data[, -c(9,10,12) ])
# Results
reg3# Coefficients only
reg4 <- lm(y~ x2)
summary(reg3)
# Using variable groups
reg3 <- lm(y ~ x)
reg4 <- lm(y~ x2)
# Results
reg3# Coefficients only
summary(reg3)
summ(reg3)
summ(reg4)
##### Confidence Intervals for Coefficients####
confint(reg3)
hist(residuals(reg3))  # Histogram of residuals
## same
hist(resid(reg3))
hist(resid(reg1))
hist(resid(reg1))
install.packages("lmtest")
library(lmtest)
bptest(reg1)
bptest(reg3)
plot (resid(reg1))
plot (resid(reg1))
# Clear environment
rm(list = ls())
# Clear console
cat("\014")  # ctrl+L
# Installs pacman ("package manager") if needed
if (!require("pacman")) install.packages("pacman")
# Use pacman to load add-on packages as desired
pacman::p_load(datasets, pacman, caret, lars, tidyverse, jtools, huxtable,
officer, flextable, texreg, AER, pastecs, openxlsx)
###### get to know the data we want to use #####
?USJudgeRatings
head(USJudgeRatings)
summary(USJudgeRatings)
Mydata1 <- USJudgeRatings %>%
rename(contacts = CONT)
# we can also select a number of variables within a dataset
Mydata2 <- USJudgeRatings %>%
select(CONT,INTG,RTEN) %>%
rename(contacts = CONT)%>%
rename(worthy= RTEN)
### Let's remove these changed datasets:
rm(Mydata1,Mydata2)
tablestat<-stat.desc(USJudgeRatings) # you need to have "pastecs" package to be able to use this function
reg1 <- lm(USJudgeRatings$RTEN ~ USJudgeRatings$CONT)
reg1 <- lm(RTEN ~ CONT, data = USJudgeRatings)
summary(reg1)
summ(reg1) # from jtools package
summ(reg1) # from jtools package
summary(reg1)
summ(reg1) # from jtools package
# visualization of the regression line
plot(RTEN ~ CONT, data = USJudgeRatings)
abline(reg1)
##############################################
# If we want to remove intercept from the linear regression model:
reg0 <- lm(RTEN ~ CONT+0, data = USJudgeRatings)
summary(reg0)
plot(RTEN ~ CONT, data = USJudgeRatings)
abline(reg0)
summary(reg2 <- lm(RTEN ~ CONT + INTG + DMNR + DILG + CFMG +
DECI + PREP + FAMI + ORAL + WRIT + PHYS,
data = USJudgeRatings))
summ(reg2)
?export_summs
export_summs(reg2)
# In order to export your regression tables in a clean format
export_summs(reg1, reg2) # to have publishable Table output
library(datasets)  # Load base packages manually
# Installs pacman ("package manager") if needed
if (!require("pacman")) install.packages("pacman")
# Use pacman to load add-on packages as desired
pacman::p_load(pacman, caret, lars, tidyverse, rio, psych,
GGally,arm,dplyr, huxtable, officer, flextable,
jtools, plyr, vtable)
########## importation des données
infidelite <- import(file.choose())
describe(infidelite)
summary(infidelite)
###### Correlation matrix
c2 <- cor(infidelite$Y, infidelite$Z3)
(Cl <- cor(infidelite)) # correlation Matrix
# You need GGally package
?ggcorr
ggcorr(infidelite, method = c("pairwise", "pearson")) #to visualize correlation
ggcorr(infidelite, method = c("pairwise", "pearson"),
nbreaks = NULL, digits = 3, low = "blue",
mid = "white", high = "red",
geom = "tile", label = T,
label_alpha = F)
#barchart (infidelite$Y,infidelite$Z1)
plot(infidelite) # not interesting
ggcorr(infidelite, method = c("pairwise", "pearson"),
nbreaks = NULL, digits = 3, low = "blue",
mid = "white", high = "red",
geom = "tile", label = T,
label_alpha = F)
#barchart (infidelite$Y,infidelite$Z1)
plot(infidelite) # not interesting
attach(infidelite)
plot(Z2,Y)  # age and infidelity, it does not look nice either
hist(infidelite$Y)
(Yfrequency <- table(Y))
barplot(Yfrequency)  # the bar chart does a better job than the histogram
summary(regA <- lm(Y ~ Z1 + Z2 + Z3 + Z4 + Z5 + Z6 + Z7 + Z8 ))
summary(regB <- lm(Y ~ Z1  + Z3 + Z4 + Z5 + Z7 + Z8 )) # to avoid multicollinearity
summary(regB <- lm(Y ~ Z1  + Z3 + Z4 + Z5 + Z7 + Z8 )) # to avoid multicollinearity
summary(regA <- lm(Y ~ Z1 + Z2 + Z3 + Z4 + Z5 + Z6 + Z7 + Z8 ))
# nice regression table:
export_summs(regA, regB)
# nice regression table:
export_summs(regA, regB)
t.test(Y ~ Z1)  # (mean(Y)for Z1=1 ) == (mean(Y)for Z1=0 )
t.test(Y ~ Z4)
detach(infidelite)
detach(infidelite)
# Clear plots
dev.off()  # But only if there IS a plot
# Clear environment
rm(list = ls())
# Clear packages
p_unload(all)  # Remove all add-ons
# Clear console
cat("\014")  # ctrl+L
# Use pacman to load add-on packages as desired
pacman::p_load(pacman, caret, lars, tidyverse, rio, psych, GGally,arm, AER)
data("Journals", package = "AER")
barchart (infidelite$Y,infidelite$Z1)
########## importation des données
infidelite <- import(file.choose())
########## importation des données
infidelite <- import(file.choose())
###### Correlation matrix
c2 <- cor(infidelite$Y, infidelite$Z3)
(Cl <- cor(infidelite)) # correlation Matrix
barchart (infidelite$Y,infidelite$Z1)
# CLEAN UP #################################################
# Clear plots
dev.off()  # But only if there IS a plot
# Clear environment
rm(list = ls())
# Clear packages
p_unload(all)  # Remove all add-ons
# Clear console
cat("\014")  # ctrl+L
# Use pacman to load add-on packages as desired
pacman::p_load(pacman, caret, lars, tidyverse, rio, psych, GGally,arm, AER)
data("Journals", package = "AER")
journals <- Journals[, c("subs", "price")]
journals$citeprice <- Journals$price/Journals$citations # creation of a new variable
summary(journals)
attach(journals)
plot(subs,citeprice)
plot(subs ~ citeprice, data = journals) # I suspect a nonlinear relationship
plot(log(subs) ~ log(citeprice), data = journals) # Using log-linearisation
?rug
rug(log(citeprice))
rug(log(subs), side = 2)
plot(log(subs) ~ log(citeprice), data = Journals, pch = 20,
col = "blue", ylim = c(0, 8), xlim = c(-7, 4),
main = "Library subscriptions")  #changing the graphical parameters
reg_jour <- lm(log(subs) ~ log(citeprice), data = journals)
?abline
abline(reg_jour) # Add straight line
summary(reg_jour)
plot(reg_jour) # plotting a regression!
plot(reg_jour, which = 1) # Plots are also accessible individually,
##  using ggplot for sophisticated graphs
?ggplot
Journals %>%
ggplot(aes(x=citations, y=price)) + geom_point()
Journals %>%
ggplot(aes(x=log(subs),
y=log(price)-log(citations),
col=society,
size=pages)) + geom_point() # A figure
Journals %>%
filter(foundingyear>=1970)%>%
ggplot(aes(x=log(subs),
y=log(price)-log(citations),
col=field, size=pages)) + geom_point()
pdf("~/Desktop/testfile.pdf", height = 5, width = 6)
#### bonus: to have an overview of the charts in R
demo("graphics")
#### bonus: to have an overview of the charts in R
demo("graphics")
#### bonus: to have an overview of the charts in R
demo("graphics")
dev.off()
#### bonus: to have an overview of the charts in R
demo("graphics")
# Clear plots
dev.off()  # But only if there IS a plot
# Clear environment
rm(list = ls())
# Clear packages
p_unload(all)  # Remove all add-ons
# Clear console
cat("\014")  # ctrl+L
pacman::p_load(pacman, caret, lars, tidyverse, rio, psych, GGally,arm, AER)
data("CPS1988", package = "AER")
dim(CPS1988)
summary(CPS1988)
describe(CPS1988)
attach(CPS1988)
ethnicityfreq <- table(ethnicity)
barchart (ethnicityfreq)
View(CPS1988)
summary(CPS1988)
summary(CPS1988)
describe(CPS1988)
attach(CPS1988)
ethnicityfreq <- table(ethnicity)
barchart (ethnicityfreq)
hist(experience)
# t-test to see whether the mean of wage is statistically different for african Americans or not
t.test(wage ~ ethnicity)
CPS1988 %>%
ggplot(aes(x=education, y=wage , col=parttime))+ geom_point()
CPS1988 %>%
filter(parttime=="no") %>%
filter(experience > 10) %>%
ggplot(aes(x=education, y=wage , col=ethnicity)) + geom_point()
hist(education)
CPS1988 %>%
filter(parttime=="no") %>%
filter(experience > 10) %>%
filter(education > 12) %>%
ggplot(aes(x=education, y=wage , col=ethnicity)) + geom_point()
detach(CPS1988)
#If no other transformation is required: I() can be used, it returns its argument “as is”
reg1 <- lm(log(wage) ~ experience + I(experience^2) +
education + ethnicity, data = CPS1988)
summary(reg1)
summary(reg1)
# Clear plots
dev.off()  # But only if there IS a plot
attach(CPS1988)
plot(wage,experience)
plot(log(wage),log(experience))
reg1 <- lm(log(wage) ~ experience + I(experience^2) +
education + ethnicity)
summary(reg1)
detach(CPS1988)
# Clear plots
dev.off()  # But only if there IS a plot
# Clear environment
rm(list = ls())
# Clear packages
p_unload(all)  # Remove all add-ons
# Clear console
cat("\014")  # ctrl+L
# Use pacman to load add-on packages as desired
pacman::p_load(pacman, caret, lars, tidyverse, rio, psych,
GGally,arm,dplyr, huxtable, officer, flextable,
jtools, plyr, margins)
#-- import the file entitled "infidelité_probit.xlsx"
infidelite <- import(file.choose())
head(infidelite)
head(infidelite)
# Define variables
Y <- cbind(Affaire)
X <- cbind(Z1, Z3, Z4, Z5, Z7, Z8)
# Descriptive statistics
summary(Y)
summary(X)
attach(infidelite)
head(infidelite)
# Define variables
Y <- cbind(Affaire)
X <- cbind(Z1, Z3, Z4, Z5, Z7, Z8)
# Descriptive statistics
summary(Y)
summary(X)
table(Y)
#-- to calculate the percentage of infidel people (25%)
table(Y)/sum(table(Y))
sum(table(Y))
# Regression coefficients
olsreg <- lm(Y ~ X)
summary(olsreg)
# Logit model coefficients
logit<- glm(Y ~ X, family=binomial (link = "logit"))
summary(logit)
# Logit model odds ratios: odds ratio is a statistic that quantifies the strength of the association between two events
exp(logit$coefficients)
# Probit model coefficients
probit<- glm(Y ~ X, family=binomial (link="probit"))
# Probit model coefficients
probit<- glm(Y ~ X, family=binomial (link="probit"))
summary(probit)
# Logit model average marginal effects  (Marginal effects can be calculated for each observation point, common calculations are to calculate the marginal effects at the mean or better--> calculate the average marginal effects)
LogitScalar <- mean(dlogis(predict(logit, type = "link")))
LogitScalar * coef(logit)
# Clear environment
rm(list = ls())
# Clear packages
p_unload(all)  # Remove all add-ons
# Clear console
cat("\014")  # ctrl+L
install.packages("WDI")
library(WDI)
# Search for GDP and population indicators
gdp_indicators <- WDIsearch("GDP")
View(gdp_indicators)
pop_indicators <- WDIsearch("population")
View(pop_indicators)
# View the search results
head(gdp_indicators)
head(pop_indicators)
View(gdp_indicators)
# Download data for GDP and population
data <- WDI(indicator = c("NY.GDP.MKTP.CD", "SP.POP.TOTL", "EG.ELC.ACCS.ZS", "SL.UEM.TOTL.ZS"),
country = c("all"),  # Specify country codes or use "all" for all countries
start = 2000,                         # Start year
end = 2022)                           # End year
# Download data for GDP and population
data <- WDI(indicator = c("NY.GDP.MKTP.CD", "SP.POP.TOTL", "EG.ELC.ACCS.ZS", "SL.UEM.TOTL.ZS"),
country = c("all"),  # Specify country codes or use "all" for all countries
start = 2000,                         # Start year
end = 2022)                           # End year
View(data)
## Lets adopt easier names for the variables
require(GGally)
## Lets adopt easier names for the variables
require(GGally)
## Lets adopt easier names for the variables
require(GGally)
require(tidyverse)
panel_data <- panel_data %>%
rename(gdp = NY.GDP.MKTP.CD) %>%
rename(pop = SP.POP.TOTL) %>%
rename(unemp = SL.UEM.TOTL.ZS) %>%
rename(electricity = EG.ELC.ACCS.ZS)
# Convert the data to a panel data frame or
# Declare panel structure:iso2c = individual, year = time
panel_data <- pdata.frame(data, index = c("iso2c", "year"))
# Check the panel data structure
head(panel_data)
# Convert the data to a panel data frame or
# Declare panel structure:iso2c = individual, year = time
panel_data <- pdata.frame(data, index = c("iso2c", "year"))
# To work with Panel data, first, Install and load the plm package if not already done
install.packages("plm")
library(plm)
# Convert the data to a panel data frame or
# Declare panel structure:iso2c = individual, year = time
panel_data <- pdata.frame(data, index = c("iso2c", "year"))
# Check the panel data structure
head(panel_data)
## Lets adopt easier names for the variables
require(GGally)
require(tidyverse)
panel_data <- panel_data %>%
rename(gdp = NY.GDP.MKTP.CD) %>%
rename(pop = SP.POP.TOTL) %>%
rename(unemp = SL.UEM.TOTL.ZS) %>%
rename(electricity = EG.ELC.ACCS.ZS)
##### Summary statistics for the entire panel data
summary(panel_data)
summary(pooled_model)
# 1. Pooled OLS regression
pooled_model <- plm(gdp ~ pop + electricity , data = panel_data, model = "pooling")
summary(pooled_model)
# 2. Fixed Effects regression (within estimator)
fe_model <- plm(gdp ~ pop + electricity, data = panel_data,  model = "within")
summary(fe_model)
summary(re_model)
# 3. Random Effects regression
# The random effects model assumes that the entity-specific effects are random and uncorrelated with the explanatory variables.
re_model <- plm(gdp ~ pop + electricity, data = panel_data,  model ="random")
summary(re_model)
View(panel_data)
# 1. Pooled OLS regression
pooled_model <- plm(gdp ~ pop + electricity + unemp, data = panel_data, model = "pooling")
summary(pooled_model)
# 1. Pooled OLS regression
pooled_model <- plm(gdp ~ pop + electricity, data = panel_data, model = "pooling")
summary(pooled_model)
# 2. Fixed Effects regression (within estimator)
fe_model <- plm(gdp ~ pop + electricity, data = panel_data,  model = "within")
summary(fe_model)
# 3. Random Effects regression
# The random effects model assumes that the entity-specific effects are random and uncorrelated with the explanatory variables.
re_model <- plm(gdp ~ pop + electricity, data = panel_data,  model ="random")
summary(re_model)
# Clear environment
rm(list = ls())
u
# Load the built-in Grunfeld panel dataset
data("Grunfeld", package = "plm")
?Grunfeld
View(Grunfeld)
# Take a look at the structure
head(Grunfeld)
# Declare panel structure: firm = individual, year = time
pdata <- pdata.frame(Grunfeld, index = c("firm", "year"))
# Declare panel structure: firm = individual, year = time
pdata <- pdata.frame(Grunfeld, index = c("firm", "year"))
# 1. Pooled OLS regression
# Regression: investment (inv) on market value (value) and capital stock (capital)
model_pooled <- plm(inv ~ value + capital, data = pdata, model = "pooling")
summary(model_pooled)
# 2. Fixed Effects regression (within estimator)
model_fe <- plm(inv ~ value + capital, data = pdata, model = "within")
summary(model_fe)
# 3. Random Effects regression
model_re <- plm(inv ~ value + capital, data = pdata, model = "random")
summary(model_re)
# 4. Hausman test to decide between Fixed vs Random Effects
phtest(model_fe, model_re)
# 4. Hausman test to decide between Fixed vs Random Effects
?phtest
phtest(model_fe, model_re)
knitr::opts_chunk$set(echo = TRUE)
pacman::p_load(datasets, psych, pastecs)
data("iris") #iris will become an object in the environment
summary(iris)
describe(iris)
# Put graphs in 3 rows and 1 column
par(mfrow = c(3, 1))
# Histograms for each species using options
hist(iris$Petal.Width [iris$Species == "setosa"],
xlim = c(0, 3),
breaks = 5,
main = "Petal Width for Setosa",
xlab = "",
col = "red")
hist(iris$Petal.Width [iris$Species == "versicolor"],
xlim = c(0, 3),
breaks = 9,
main = "Petal Width for Versicolor",
xlab = "",
col = "purple")
hist(iris$Petal.Width [iris$Species == "virginica"],
xlim = c(0, 3),
breaks = 9,
main = "Petal Width for Virginica",
xlab = "",
col = "blue")
# Restore graphic parameter
par(mfrow=c(1, 1))
library(tinytex)
pacman::p_load(datasets, psych, pastecs)
data("iris") #iris will become an object in the environment
summary(iris)
describe(iris)
# Put graphs in 3 rows and 1 column
par(mfrow = c(3, 1))
# Histograms for each species using options
hist(iris$Petal.Width [iris$Species == "setosa"],
xlim = c(0, 3),
breaks = 5,
main = "Petal Width for Setosa",
xlab = "",
col = "red")
hist(iris$Petal.Width [iris$Species == "versicolor"],
xlim = c(0, 3),
breaks = 9,
main = "Petal Width for Versicolor",
xlab = "",
col = "purple")
hist(iris$Petal.Width [iris$Species == "virginica"],
xlim = c(0, 3),
breaks = 9,
main = "Petal Width for Virginica",
xlab = "",
col = "blue")
# Restore graphic parameter
par(mfrow=c(1, 1))
(setosa_versicolor_ttest <- t.test(
iris$Sepal.Width[iris$Species == "setosa"],
iris$Sepal.Width[iris$Species == "versicolor"]
))
plot(iris$Petal.Length, iris$Petal.Width,
col = "#75197F",  # Hex code for datalab.cc red
pch = 6,         # Use solid circles for points
main = "Iris: Petal Length vs. Petal Width",
xlab = "Petal Length",
ylab = "Petal Width")
ggcorr(iris, method = c("pairwise", "pearson"),
nbreaks = NULL, digits = 3, low = "blue",
mid = "white", high = "red",
geom = "tile", label = T,
label_alpha = F)
# You need GGally package
?ggcorr
